How To: Orthorectify EMIT Data¶


EMIT Data is provided in non-orthorectified format to reduce data size. The location group contains latitude and longitude values of each pixel as well as a geometric lookup table (GLT) that can be used to orthocorrect the imagery. The GLT is an array that provides relative downtrack and crosstrack reference locations from the raw scene to facilitate fast projection of the dataset.

This notebook will walk through using the emit_tools module to do the orthorectification process as well as a in depth walkthrough of the process.

If you are planning to convert the EMIT netCDF4 files to .envi format, the reformat.py available in the emit-sds/emit-utils repository can be used to orthorectify the imagery during the conversion process.

Requirements:

  • A NASA Earthdata Login account is required to download EMIT data
  • Selected the emit_tutorials environment as the kernel for this notebook.
    • For instructions on setting up the environment, follow the instructions in the Prerequisites section of the README.md included in the repository.
  • Downloaded the necessary EMIT files to the ../data/ folder.
    • For download instructions see the Prerequisites section of the README.md included in the repository.

Learning Objectives

  • How to open an EMIT file as an xarray.Dataset
  • How to orthorectify an EMIT file

Import the required Python libraries.

In [1]:
# Import Packages
import os
import netCDF4 as nc
from osgeo import gdal
import numpy as np
import xarray as xr
import hvplot.xarray
import holoviews as hv
import sys
sys.path.append('../modules/')
from emit_tools import emit_xarray

Next, set the path to an EMIT dataset as an object. In this example we use an L2A Reflectance file.

In [2]:
fp = '../data/EMIT_L2A_RFL_001_20220903T163129_2224611_012.nc'

Import the emit_xarray function from the emit_tools module and use it to read in the data.

The emit_xarray function will open and and place data from the groups of an EMIT dataset into an xarray.Dataset. It includes an ortho option which is set to True by default. This function reads the root group, location group, and sensor_band_parameters groups from the EMIT dataset. It then flattens all of the variables contained within those 3 groups into a single xarray_dataset. When ortho=True, it uses the GLT layers to build a lat/lon grid and places the values from the reflectance root group into the grid. The lat and lon grid is used as dimensional coordinates in the output xarray.dataset while the wavelengths and fwhm from the sensor_band_parameters group are used as non-dimensional coordinates. The reflectance data from the root group of the netCDF is then orthorectified and added as a variable in the xarray.Dataset. This function also includes an optional quality_mask parameter, which can be used to mask poor quality data. For more on masking, see How to use EMIT Quality Data.

Read in a dataset using the emit_xarray function.

In [3]:
ds = emit_xarray(fp)
ds
Out[3]:
<xarray.Dataset>
Dimensions:      (latitude: 2009, longitude: 2353, bands: 285)
Coordinates:
  * latitude     (latitude) float64 -39.31 -39.31 -39.31 ... -40.39 -40.4 -40.4
  * longitude    (longitude) float64 -62.51 -62.51 -62.51 ... -61.24 -61.24
    wavelengths  (bands) float32 381.0 388.4 395.8 ... 2.486e+03 2.493e+03
    fwhm         (bands) float32 8.415 8.415 8.415 8.415 ... 8.806 8.807 8.809
    spatial_ref  int32 0
Dimensions without coordinates: bands
Data variables:
    reflectance  (latitude, longitude, bands) float32 nan nan nan ... nan nan
Attributes: (12/38)
    ncei_template_version:             NCEI_NetCDF_Swath_Template_v2.0
    summary:                           The Earth Surface Mineral Dust Source ...
    keywords:                          Imaging Spectroscopy, minerals, EMIT, ...
    Conventions:                       CF-1.63
    sensor:                            EMIT (Earth Surface Mineral Dust Sourc...
    instrument:                        EMIT
    ...                                ...
    southernmost_latitude:             -40.39610428069674
    spatialResolution:                 0.000542232520256367
    spatial_ref:                       GEOGCS["WGS 84",DATUM["WGS_1984",SPHER...
    geotransform:                      [-6.25120945e+01  5.42232520e-04 -0.00...
    day_night_flag:                    Day
    title:                             EMIT L2A Surface Reflectance 60 m V001
xarray.Dataset
    • latitude: 2009
    • longitude: 2353
    • bands: 285
    • latitude
      (latitude)
      float64
      -39.31 -39.31 ... -40.4 -40.4
      long_name :
      Latitude (WGS-84)
      units :
      degrees north
      array([-39.306759, -39.307301, -39.307844, ..., -40.394478, -40.39502 ,
             -40.395562])
    • longitude
      (longitude)
      float64
      -62.51 -62.51 ... -61.24 -61.24
      long_name :
      Longitude (WGS-84)
      units :
      degrees east
      array([-62.512095, -62.511552, -62.51101 , ..., -61.237848, -61.237306,
             -61.236764])
    • wavelengths
      (bands)
      float32
      381.0 388.4 ... 2.486e+03 2.493e+03
      long_name :
      Wavelength Centers
      units :
      nm
      array([ 381.00558,  388.4092 ,  395.81583,  403.2254 ,  410.638  ,
              418.0536 ,  425.47214,  432.8927 ,  440.31726,  447.7428 ,
              455.17035,  462.59888,  470.0304 ,  477.46292,  484.89743,
              492.33292,  499.77142,  507.2099 ,  514.6504 ,  522.0909 ,
              529.5333 ,  536.9768 ,  544.42126,  551.8667 ,  559.3142 ,
              566.7616 ,  574.20905,  581.6585 ,  589.108  ,  596.55835,
              604.0098 ,  611.4622 ,  618.9146 ,  626.36804,  633.8215 ,
              641.2759 ,  648.7303 ,  656.1857 ,  663.6411 ,  671.09753,
              678.5539 ,  686.0103 ,  693.4677 ,  700.9251 ,  708.38354,
              715.84094,  723.2993 ,  730.7587 ,  738.2171 ,  745.6765 ,
              753.1359 ,  760.5963 ,  768.0557 ,  775.5161 ,  782.97754,
              790.4379 ,  797.89935,  805.36176,  812.8232 ,  820.2846 ,
              827.746  ,  835.2074 ,  842.66986,  850.1313 ,  857.5937 ,
              865.0551 ,  872.5176 ,  879.98004,  887.44147,  894.90393,
              902.3664 ,  909.82886,  917.2913 ,  924.7538 ,  932.21625,
              939.6788 ,  947.14026,  954.6027 ,  962.0643 ,  969.5268 ,
              976.9883 ,  984.4498 ,  991.9114 ,  999.37286, 1006.8344 ,
             1014.295  , 1021.7566 , 1029.2172 , 1036.6777 , 1044.1383 ,
             1051.5989 , 1059.0596 , 1066.5201 , 1073.9797 , 1081.4404 ,
             1088.9    , 1096.3597 , 1103.8184 , 1111.2781 , 1118.7368 ,
      ...
             1796.4385 , 1803.8701 , 1811.3008 , 1818.7314 , 1826.1611 ,
             1833.591  , 1841.0206 , 1848.4495 , 1855.8773 , 1863.3052 ,
             1870.733  , 1878.16   , 1885.5869 , 1893.013  , 1900.439  ,
             1907.864  , 1915.2892 , 1922.7133 , 1930.1375 , 1937.5607 ,
             1944.9839 , 1952.4071 , 1959.8295 , 1967.2518 , 1974.6732 ,
             1982.0946 , 1989.515  , 1996.9355 , 2004.355  , 2011.7745 ,
             2019.1931 , 2026.6118 , 2034.0304 , 2041.4471 , 2048.865  ,
             2056.2808 , 2063.6965 , 2071.1123 , 2078.5273 , 2085.9421 ,
             2093.3562 , 2100.769  , 2108.1821 , 2115.5942 , 2123.0063 ,
             2130.4175 , 2137.8289 , 2145.239  , 2152.6482 , 2160.0576 ,
             2167.467  , 2174.8755 , 2182.283  , 2189.6904 , 2197.097  ,
             2204.5034 , 2211.9092 , 2219.3147 , 2226.7195 , 2234.1233 ,
             2241.5269 , 2248.9297 , 2256.3328 , 2263.7346 , 2271.1365 ,
             2278.5376 , 2285.9387 , 2293.3386 , 2300.7378 , 2308.136  ,
             2315.5342 , 2322.9326 , 2330.3298 , 2337.7263 , 2345.1216 ,
             2352.517  , 2359.9126 , 2367.3071 , 2374.7007 , 2382.0935 ,
             2389.486  , 2396.878  , 2404.2695 , 2411.6604 , 2419.0513 ,
             2426.4402 , 2433.8303 , 2441.2183 , 2448.6064 , 2455.9944 ,
             2463.3816 , 2470.7678 , 2478.153  , 2485.5386 , 2492.9238 ],
            dtype=float32)
    • fwhm
      (bands)
      float32
      8.415 8.415 8.415 ... 8.807 8.809
      long_name :
      Full Width at Half Max
      units :
      nm
      array([8.415, 8.415, 8.415, 8.415, 8.417, 8.418, 8.419, 8.421, 8.422,
             8.424, 8.425, 8.426, 8.428, 8.429, 8.431, 8.432, 8.433, 8.435,
             8.436, 8.438, 8.439, 8.44 , 8.442, 8.443, 8.445, 8.446, 8.447,
             8.449, 8.45 , 8.452, 8.453, 8.454, 8.456, 8.457, 8.459, 8.46 ,
             8.461, 8.463, 8.464, 8.466, 8.467, 8.468, 8.47 , 8.471, 8.473,
             8.474, 8.475, 8.477, 8.478, 8.48 , 8.481, 8.482, 8.484, 8.485,
             8.487, 8.488, 8.489, 8.491, 8.492, 8.494, 8.495, 8.496, 8.498,
             8.499, 8.501, 8.502, 8.503, 8.505, 8.506, 8.508, 8.509, 8.51 ,
             8.512, 8.513, 8.515, 8.516, 8.517, 8.519, 8.52 , 8.522, 8.523,
             8.524, 8.526, 8.527, 8.529, 8.53 , 8.531, 8.533, 8.534, 8.536,
             8.537, 8.538, 8.54 , 8.541, 8.543, 8.544, 8.545, 8.547, 8.548,
             8.55 , 8.551, 8.552, 8.554, 8.555, 8.557, 8.558, 8.559, 8.561,
             8.562, 8.564, 8.565, 8.566, 8.568, 8.569, 8.571, 8.572, 8.573,
             8.575, 8.576, 8.578, 8.579, 8.58 , 8.582, 8.583, 8.585, 8.586,
             8.587, 8.589, 8.59 , 8.592, 8.593, 8.594, 8.596, 8.597, 8.599,
             8.6  , 8.601, 8.603, 8.604, 8.606, 8.607, 8.608, 8.61 , 8.611,
             8.613, 8.614, 8.615, 8.617, 8.618, 8.62 , 8.621, 8.622, 8.624,
             8.625, 8.627, 8.628, 8.629, 8.631, 8.632, 8.634, 8.635, 8.636,
             8.638, 8.639, 8.641, 8.642, 8.643, 8.645, 8.646, 8.648, 8.649,
             8.65 , 8.652, 8.653, 8.655, 8.656, 8.657, 8.659, 8.66 , 8.662,
             8.663, 8.664, 8.666, 8.667, 8.669, 8.67 , 8.671, 8.673, 8.674,
             8.676, 8.677, 8.678, 8.68 , 8.681, 8.683, 8.684, 8.685, 8.687,
             8.688, 8.69 , 8.691, 8.692, 8.694, 8.695, 8.697, 8.698, 8.699,
             8.701, 8.702, 8.704, 8.705, 8.706, 8.708, 8.709, 8.711, 8.712,
             8.714, 8.715, 8.716, 8.718, 8.719, 8.721, 8.722, 8.723, 8.725,
             8.726, 8.727, 8.729, 8.73 , 8.732, 8.733, 8.734, 8.736, 8.737,
             8.739, 8.74 , 8.741, 8.743, 8.744, 8.746, 8.747, 8.748, 8.75 ,
             8.751, 8.753, 8.754, 8.755, 8.757, 8.758, 8.76 , 8.761, 8.763,
             8.764, 8.765, 8.767, 8.768, 8.77 , 8.771, 8.772, 8.774, 8.775,
             8.777, 8.778, 8.779, 8.781, 8.782, 8.783, 8.785, 8.786, 8.788,
             8.789, 8.79 , 8.792, 8.793, 8.795, 8.796, 8.797, 8.799, 8.8  ,
             8.802, 8.803, 8.804, 8.806, 8.807, 8.809], dtype=float32)
    • spatial_ref
      ()
      int32
      0
      crs_wkt :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      grid_mapping_name :
      latitude_longitude
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      array(0)
    • reflectance
      (latitude, longitude, bands)
      float32
      nan nan nan nan ... nan nan nan nan
      long_name :
      Surface Reflectance
      units :
      unitless
      array([[[nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              ...,
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan]],
      
             [[nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              ...,
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan]],
      
             [[nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              ...,
      ...
              ...,
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan]],
      
             [[nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              ...,
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan]],
      
             [[nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              ...,
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan],
              [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(Float64Index([  -39.3067591475017, -39.307301380021954, -39.307843612542214,
                     -39.30838584506247,  -39.30892807758273,  -39.30947031010298,
                     -39.31001254262324, -39.310554775143494, -39.311097007663754,
                     -39.31163924018401,
                    ...
                     -40.39068195549418, -40.391224188014434, -40.391766420534694,
                     -40.39230865305495,  -40.39285088557521,  -40.39339311809546,
                     -40.39393535061572,  -40.39447758313597,  -40.39501981565623,
                     -40.39556204817649],
                   dtype='float64', name='latitude', length=2009))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([  -62.5120945327963, -62.511552300276044, -62.511010067755784,
                     -62.51046783523553,  -62.50992560271527,  -62.50938337019502,
                     -62.50884113767476, -62.508298905154504, -62.507756672634244,
                     -62.50721444011399,
                    ...
                     -61.24164373783563,  -61.24110150531537,  -61.24055927279512,
                     -61.24001704027486, -61.239474807754604, -61.238932575234344,
                     -61.23839034271409,  -61.23784811019384,  -61.23730587767358,
                    -61.236763645153324],
                   dtype='float64', name='longitude', length=2353))
  • ncei_template_version :
    NCEI_NetCDF_Swath_Template_v2.0
    summary :
    The Earth Surface Mineral Dust Source Investigation (EMIT) is an Earth Ventures-Instrument (EVI-4) Mission that maps the surface mineralogy of arid dust source regions via imaging spectroscopy in the visible and short-wave infrared (VSWIR). Installed on the International Space Station (ISS), the EMIT instrument is a Dyson imaging spectrometer that uses contiguous spectroscopic measurements from 410 to 2450 nm to resolve absoprtion features of iron oxides, clays, sulfates, carbonates, and other dust-forming minerals. During its one-year mission, EMIT will observe the sunlit Earth's dust source regions that occur within +/-52° latitude and produce maps of the source regions that can be used to improve forecasts of the role of mineral dust in the radiative forcing (warming or cooling) of the atmosphere.\n\nThis file contains L2A estimated surface reflectances and geolocation data. Reflectance estimates are created using an Optimal Estimation technique - see ATBD for details. Reflectance values are reported as fractions (relative to 1).
    keywords :
    Imaging Spectroscopy, minerals, EMIT, dust, radiative forcing
    Conventions :
    CF-1.63
    sensor :
    EMIT (Earth Surface Mineral Dust Source Investigation)
    instrument :
    EMIT
    platform :
    ISS
    processing_version :
    V1.0
    institution :
    NASA Jet Propulsion Laboratory/California Institute of Technology
    license :
    https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/
    naming_authority :
    LPDAAC
    date_created :
    2022-11-14T09:50:54Z
    keywords_vocabulary :
    NASA Global Change Master Directory (GCMD) Science Keywords
    stdname_vocabulary :
    NetCDF Climate and Forecast (CF) Metadata Convention
    creator_name :
    Jet Propulsion Laboratory/California Institute of Technology
    creator_email :
    sarah.r.lundeen@jpl.nasa.gov
    creator_url :
    https://earth.jpl.nasa.gov/emit/
    project :
    Earth Surface Mineral Dust Source Investigation
    project_url :
    https://emit.jpl.nasa.gov/
    publisher_name :
    USGS LPDAAC
    publisher_url :
    https://lpdaac.usgs.gov
    publisher_email :
    lpdaac@usgs.gov
    identifier_product_doi_authority :
    https://doi.org
    flight_line :
    emit20220903t163129_o24611_s000
    time_coverage_start :
    2022-09-03T16:31:29+0000
    time_coverage_end :
    2022-09-03T16:31:41+0000
    software_build_version :
    010603
    product_version :
    01
    history :
    PGE Input files: radiance_file=/beegfs/store/emit/ops/data/acquisitions/20220903/emit20220903t163129/l1b/emit20220903t163129_o24611_s000_l1b_rdn_b0106_v01.img, pixel_locations_file=/beegfs/store/emit/ops/data/acquisitions/20220903/emit20220903t163129/l1b/emit20220903t163129_o24611_s000_l1b_loc_b0106_v01.img, observation_parameters_file=/beegfs/store/emit/ops/data/acquisitions/20220903/emit20220903t163129/l1b/emit20220903t163129_o24611_s000_l1b_obs_b0106_v01.img, surface_model_config=/beegfs/store/emit/ops/repos/emit-sds-l2a/surface/surface_20221020.json
    easternmost_longitude :
    -62.5120945327963
    northernmost_latitude :
    -39.3067591475017
    westernmost_longitude :
    -61.236221412633064
    southernmost_latitude :
    -40.39610428069674
    spatialResolution :
    0.000542232520256367
    spatial_ref :
    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
    geotransform :
    [-6.25120945e+01 5.42232520e-04 -0.00000000e+00 -3.93067591e+01 -0.00000000e+00 -5.42232520e-04]
    day_night_flag :
    Day
    title :
    EMIT L2A Surface Reflectance 60 m V001

Now plot the results of a red band.

In [4]:
# Get red band and plot example
b650 = np.nanargmin(abs(ds['wavelengths'].data-650))
ds.isel(bands=b650).hvplot.image(cmap='viridis', aspect = 'equal', frame_width=500, rasterize=True)
Out[4]:

Now open a dataset with ortho=False. The main difference here outside of the data not being orthocorrected is that there are no coordinates associated with this xarray.dataset, all data from the netCDF are included, but are treated as variables.

In [5]:
ds = emit_xarray(fp,ortho=False)
ds
Out[5]:
<xarray.Dataset>
Dimensions:      (downtrack: 1280, crosstrack: 1242, bands: 285, ortho_y: 2009,
                  ortho_x: 2353)
Dimensions without coordinates: downtrack, crosstrack, bands, ortho_y, ortho_x
Data variables:
    reflectance  (downtrack, crosstrack, bands) float32 ...
    lat          (downtrack, crosstrack) float64 ...
    lon          (downtrack, crosstrack) float64 ...
    elev         (downtrack, crosstrack) float64 ...
    glt_x        (ortho_y, ortho_x) float64 ...
    glt_y        (ortho_y, ortho_x) float64 ...
    wavelengths  (bands) float32 ...
    fwhm         (bands) float32 ...
Attributes: (12/38)
    ncei_template_version:             NCEI_NetCDF_Swath_Template_v2.0
    summary:                           The Earth Surface Mineral Dust Source ...
    keywords:                          Imaging Spectroscopy, minerals, EMIT, ...
    Conventions:                       CF-1.63
    sensor:                            EMIT (Earth Surface Mineral Dust Sourc...
    instrument:                        EMIT
    ...                                ...
    southernmost_latitude:             -40.39610428069674
    spatialResolution:                 0.000542232520256367
    spatial_ref:                       GEOGCS["WGS 84",DATUM["WGS_1984",SPHER...
    geotransform:                      [-6.25120945e+01  5.42232520e-04 -0.00...
    day_night_flag:                    Day
    title:                             EMIT L2A Surface Reflectance 60 m V001
xarray.Dataset
    • downtrack: 1280
    • crosstrack: 1242
    • bands: 285
    • ortho_y: 2009
    • ortho_x: 2353
      • reflectance
        (downtrack, crosstrack, bands)
        float32
        ...
        long_name :
        Surface Reflectance
        units :
        unitless
        [453081600 values with dtype=float32]
      • lat
        (downtrack, crosstrack)
        float64
        ...
        long_name :
        Longitude (WGS-84)
        units :
        degrees east
        [1589760 values with dtype=float64]
      • lon
        (downtrack, crosstrack)
        float64
        ...
        long_name :
        Latitude (WGS-84)
        units :
        degrees north
        [1589760 values with dtype=float64]
      • elev
        (downtrack, crosstrack)
        float64
        ...
        long_name :
        Surface Elevation
        units :
        m
        [1589760 values with dtype=float64]
      • glt_x
        (ortho_y, ortho_x)
        float64
        ...
        long_name :
        GLT Sample Lookup
        units :
        pixel location
        [4727177 values with dtype=float64]
      • glt_y
        (ortho_y, ortho_x)
        float64
        ...
        long_name :
        GLT Line Lookup
        units :
        pixel location
        [4727177 values with dtype=float64]
      • wavelengths
        (bands)
        float32
        ...
        long_name :
        Wavelength Centers
        units :
        nm
        [285 values with dtype=float32]
      • fwhm
        (bands)
        float32
        ...
        long_name :
        Full Width at Half Max
        units :
        nm
        [285 values with dtype=float32]
      • ncei_template_version :
        NCEI_NetCDF_Swath_Template_v2.0
        summary :
        The Earth Surface Mineral Dust Source Investigation (EMIT) is an Earth Ventures-Instrument (EVI-4) Mission that maps the surface mineralogy of arid dust source regions via imaging spectroscopy in the visible and short-wave infrared (VSWIR). Installed on the International Space Station (ISS), the EMIT instrument is a Dyson imaging spectrometer that uses contiguous spectroscopic measurements from 410 to 2450 nm to resolve absoprtion features of iron oxides, clays, sulfates, carbonates, and other dust-forming minerals. During its one-year mission, EMIT will observe the sunlit Earth's dust source regions that occur within +/-52° latitude and produce maps of the source regions that can be used to improve forecasts of the role of mineral dust in the radiative forcing (warming or cooling) of the atmosphere.\n\nThis file contains L2A estimated surface reflectances and geolocation data. Reflectance estimates are created using an Optimal Estimation technique - see ATBD for details. Reflectance values are reported as fractions (relative to 1).
        keywords :
        Imaging Spectroscopy, minerals, EMIT, dust, radiative forcing
        Conventions :
        CF-1.63
        sensor :
        EMIT (Earth Surface Mineral Dust Source Investigation)
        instrument :
        EMIT
        platform :
        ISS
        processing_version :
        V1.0
        institution :
        NASA Jet Propulsion Laboratory/California Institute of Technology
        license :
        https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/
        naming_authority :
        LPDAAC
        date_created :
        2022-11-14T09:50:54Z
        keywords_vocabulary :
        NASA Global Change Master Directory (GCMD) Science Keywords
        stdname_vocabulary :
        NetCDF Climate and Forecast (CF) Metadata Convention
        creator_name :
        Jet Propulsion Laboratory/California Institute of Technology
        creator_email :
        sarah.r.lundeen@jpl.nasa.gov
        creator_url :
        https://earth.jpl.nasa.gov/emit/
        project :
        Earth Surface Mineral Dust Source Investigation
        project_url :
        https://emit.jpl.nasa.gov/
        publisher_name :
        USGS LPDAAC
        publisher_url :
        https://lpdaac.usgs.gov
        publisher_email :
        lpdaac@usgs.gov
        identifier_product_doi_authority :
        https://doi.org
        flight_line :
        emit20220903t163129_o24611_s000
        time_coverage_start :
        2022-09-03T16:31:29+0000
        time_coverage_end :
        2022-09-03T16:31:41+0000
        software_build_version :
        010603
        product_version :
        01
        history :
        PGE Input files: radiance_file=/beegfs/store/emit/ops/data/acquisitions/20220903/emit20220903t163129/l1b/emit20220903t163129_o24611_s000_l1b_rdn_b0106_v01.img, pixel_locations_file=/beegfs/store/emit/ops/data/acquisitions/20220903/emit20220903t163129/l1b/emit20220903t163129_o24611_s000_l1b_loc_b0106_v01.img, observation_parameters_file=/beegfs/store/emit/ops/data/acquisitions/20220903/emit20220903t163129/l1b/emit20220903t163129_o24611_s000_l1b_obs_b0106_v01.img, surface_model_config=/beegfs/store/emit/ops/repos/emit-sds-l2a/surface/surface_20221020.json
        easternmost_longitude :
        -62.5120945327963
        northernmost_latitude :
        -39.3067591475017
        westernmost_longitude :
        -61.236221412633064
        southernmost_latitude :
        -40.39610428069674
        spatialResolution :
        0.000542232520256367
        spatial_ref :
        GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
        geotransform :
        [-6.25120945e+01 5.42232520e-04 -0.00000000e+00 -3.93067591e+01 -0.00000000e+00 -5.42232520e-04]
        day_night_flag :
        Day
        title :
        EMIT L2A Surface Reflectance 60 m V001

      We can see the difference in structure and the plotted data if the dataset is orthorectified or not.

      In [6]:
      # Get red band and plot example
      b650 = np.nanargmin(abs(ds['wavelengths'].data-650))
      ds.isel(bands=b650).hvplot.image(cmap='viridis', aspect = 'equal', frame_width=500, rasterize=True)
      
      Out[6]:

      The orthorectified dataset can also be saved as a netCDF4 output that can be reopened using the xarray.open_dataset function if so desired using the cell below.

      In [ ]:
      ds.to_netcdf('../data/example_emit_ortho_out.nc')
      # Example for Opening 
      # ds = xr.open_dataset('../data/example_emit_ortho_out.nc')
      

      2. Detailed Explanation of Orthorectification and xarray.Dataset Construction¶

      As mentioned, the orthorectification process for EMIT data utilizes a premade geometric lookup table (GLT) that is included in the location group of the file. glt_x and glt_y contain the pixel indices from the raw dataset to construct an orthocorrected array. These indices data can be pulled from the raw dataset to populate an orthocorrected array using numpy broadcasting.

      To do this, first read in the location group data.

      In [8]:
      loc = xr.open_dataset(fp, engine = 'h5netcdf', group='location')
      

      Next, build a numpy array using the GLT arrays from the location group. Use the stack function to stack the x and y GLT arrays and use the function nan_to_num to set the NaN values to the GLT_NODATA_VALUE, which is 0 for EMIT datasets. After doing this, you can check the shape of the array to see the expected height and width in pixels of the image after orthorectification.

      In [9]:
      GLT_NODATA_VALUE=0
      glt_array = np.nan_to_num(np.stack([loc['glt_x'].data,loc['glt_y'].data],axis=-1),nan=GLT_NODATA_VALUE).astype(int)
      glt_array.shape
      
      Out[9]:
      (2009, 2353, 2)

      After we have a GLT array, we need to open and create a numpy array from the reflectance dataset. Open the dataset using xarray and assign the data to an array. Here you can check the shape of the array to see the dimensions of the uncorrected data.

      In [10]:
      ds = xr.open_dataset(fp,engine = 'h5netcdf')
      ds_array = ds['reflectance'].data
      ds_array.shape
      
      Out[10]:
      (1280, 1242, 285)

      Now that we have both arrays, we need to assign a fill value for positions that have no data. We can then construct an empty array with the dimensions we desire (GLT shape) and populate it with values from the dataset. To do this we first use np.zeros to create an array of all zeros then add the fill_value to the array of zeros to change the zeros to the fill_value.

      Next we can build an array of valid locations by omitting the portions of the array containing the GLT_NODATA_VALUES.

      After that, change the indexing of the GLT array which is one based to zero based by subtracting 1 from them.

      Lastly, we can use broadcasting/indexing to pull through the values we want from the dataset into our new output array.

      In [11]:
      # Build Output Dataset
      fill_value = -9999
      out_ds = np.zeros((glt_array.shape[0], glt_array.shape[1], ds_array.shape[-1]), dtype=np.float32) + fill_value
      valid_glt = np.all(glt_array != GLT_NODATA_VALUE, axis=-1)
      # Adjust for One based Index
      glt_array[valid_glt] -= 1 
      # Use indexing/broadcasting to populate array cells with 0 values
      out_ds[valid_glt, :] = ds_array[glt_array[valid_glt, 1], glt_array[valid_glt, 0], :]
      out_ds.shape
      
      Out[11]:
      (2009, 2353, 285)

      After this step, we have an array of values that is orthorectified, but if we want to have a grid of lat/lon values we still need to calculate them using the geotransform. We can retrieve this information from the dataset attributes (ds.attrs['geotransform']).

      In [12]:
      GT = ds.geotransform
      GT
      
      Out[12]:
      array([-6.25120945e+01,  5.42232520e-04, -0.00000000e+00, -3.93067591e+01,
             -0.00000000e+00, -5.42232520e-04])

      The geotransform consists of 2 formulas with 6 coefficients used to calculate the xy-grid to project the dataset.

      x_geo = GT[0] + x * GT[1] + y * GT[2]

      y_geo = GT[3] + x * GT[4] + y * GT[5]

      GT[0] - The x-coordinate of the upper-left corner of the upper-left pixel
      GT[1] - W-E pixel width
      GT[2] - Row rotation (zero for North up images)
      GT[3] - The y-coordinate of the upper-left corner of the upper-left pixel
      GT[4] - Column rotation (zero for North up images)
      GT[5] - N-S pixel height (negative value for a north-up image)

      Create empty arrays for the x and y (lon and lat) based on the dimensions of the dataset.

      In [13]:
      # Create Array for Lat and Lon and fill
      dim_x = loc.glt_x.shape[1]
      dim_y = loc.glt_x.shape[0]
      lon = np.zeros(dim_x)
      lat = np.zeros(dim_y)
      

      After we have an empty array with the correct dimensions, we can populate it with values using the geotransform formula to build a lat/lon grid for plotting the orthorectified data. We can remove the GT[2] and GT[4] terms from the equation since our orthorectified image is North up. Calclutate the latitude and longitude for each row (x_dim) and column (y_dim) of the dataset.

      In [14]:
      for x in np.arange(dim_x):
          x_geo = GT[0] + x * GT[1]
          lon[x] = x_geo
      for y in np.arange(dim_y):
          y_geo = GT[3] + y * GT[5]
          lat[y] = y_geo
      
      lat, lon
      
      Out[14]:
      (array([-39.30675915, -39.30730138, -39.30784361, ..., -40.39447758,
              -40.39501982, -40.39556205]),
       array([-62.51209453, -62.5115523 , -62.51101007, ..., -61.23784811,
              -61.23730588, -61.23676365]))

      As a last step we can place this and remaining information from the EMIT dataset, such as attributes and wavelengths into an xarray.Dataset if desired. To do this, build dictionaries for the coordinates and variables components of the dataset. In this example, the variable would be reflectance. First, we'll want to read in the sensor_band_parameters group to get the wavelength and fwhm info, then build the dictionaries.

      In [15]:
      wvl = xr.open_dataset(fp, engine = 'h5netcdf', group='sensor_band_parameters')
      coords = {'lat':(['lat'],lat), 'lon':(['lon'],lon), **wvl.variables} ## ** upacks the existing dictionary from the wvl dataset.
      data_vars = {'reflectance':(['lat','lon','bands'], out_ds)}
      

      Now build an xarray.Dataset using the dictionaries, then add attributes to the dataset from the non-orthocorrected data. We also can add the crs using the rio.write_crs function, to allow rasterio and rioxarray to recognize the CRS for future use.

      Note: the wavelength and fwhm are included as non-dimensional coordinates because they are independent of Lat/Lon.

      In [16]:
      out_xr = xr.Dataset(data_vars=data_vars, coords=coords, attrs= ds.attrs)
      out_xr['reflectance'].attrs = ds['reflectance'].attrs
      out_xr.coords['lat'].attrs = loc['lat'].attrs
      out_xr.coords['lon'].attrs = loc['lon'].attrs
      out_xr.rio.write_crs(ds.spatial_ref,inplace=True) # Add CRS in easily recognizable format
      out_xr
      
      Out[16]:
      <xarray.Dataset>
      Dimensions:      (lat: 2009, lon: 2353, bands: 285)
      Coordinates:
        * lat          (lat) float64 -39.31 -39.31 -39.31 ... -40.39 -40.4 -40.4
        * lon          (lon) float64 -62.51 -62.51 -62.51 ... -61.24 -61.24 -61.24
          wavelengths  (bands) float32 ...
          fwhm         (bands) float32 ...
          spatial_ref  int32 0
      Dimensions without coordinates: bands
      Data variables:
          reflectance  (lat, lon, bands) float32 -9.999e+03 -9.999e+03 ... -9.999e+03
      Attributes: (12/38)
          ncei_template_version:             NCEI_NetCDF_Swath_Template_v2.0
          summary:                           The Earth Surface Mineral Dust Source ...
          keywords:                          Imaging Spectroscopy, minerals, EMIT, ...
          Conventions:                       CF-1.63
          sensor:                            EMIT (Earth Surface Mineral Dust Sourc...
          instrument:                        EMIT
          ...                                ...
          southernmost_latitude:             -40.39610428069674
          spatialResolution:                 0.000542232520256367
          spatial_ref:                       GEOGCS["WGS 84",DATUM["WGS_1984",SPHER...
          geotransform:                      [-6.25120945e+01  5.42232520e-04 -0.00...
          day_night_flag:                    Day
          title:                             EMIT L2A Surface Reflectance 60 m V001
      xarray.Dataset
        • lat: 2009
        • lon: 2353
        • bands: 285
        • lat
          (lat)
          float64
          -39.31 -39.31 ... -40.4 -40.4
          long_name :
          Longitude (WGS-84)
          units :
          degrees east
          array([-39.306759, -39.307301, -39.307844, ..., -40.394478, -40.39502 ,
                 -40.395562])
        • lon
          (lon)
          float64
          -62.51 -62.51 ... -61.24 -61.24
          long_name :
          Latitude (WGS-84)
          units :
          degrees north
          array([-62.512095, -62.511552, -62.51101 , ..., -61.237848, -61.237306,
                 -61.236764])
        • wavelengths
          (bands)
          float32
          ...
          long_name :
          Wavelength Centers
          units :
          nm
          [285 values with dtype=float32]
        • fwhm
          (bands)
          float32
          ...
          long_name :
          Full Width at Half Max
          units :
          nm
          [285 values with dtype=float32]
        • spatial_ref
          ()
          int32
          0
          crs_wkt :
          GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
          semi_major_axis :
          6378137.0
          semi_minor_axis :
          6356752.314245179
          inverse_flattening :
          298.257223563
          reference_ellipsoid_name :
          WGS 84
          longitude_of_prime_meridian :
          0.0
          prime_meridian_name :
          Greenwich
          geographic_crs_name :
          WGS 84
          grid_mapping_name :
          latitude_longitude
          spatial_ref :
          GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
          array(0)
        • reflectance
          (lat, lon, bands)
          float32
          -9.999e+03 ... -9.999e+03
          long_name :
          Surface Reflectance
          units :
          unitless
          array([[[-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  ...,
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.]],
          
                 [[-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  ...,
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.]],
          
                 [[-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  ...,
          ...
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.]],
          
                 [[-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  ...,
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.]],
          
                 [[-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  ...,
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.],
                  [-9999., -9999., -9999., ..., -9999., -9999., -9999.]]],
                dtype=float32)
        • lat
          PandasIndex
          PandasIndex(Float64Index([  -39.3067591475017, -39.307301380021954, -39.307843612542214,
                         -39.30838584506247,  -39.30892807758273,  -39.30947031010298,
                         -39.31001254262324, -39.310554775143494, -39.311097007663754,
                         -39.31163924018401,
                        ...
                         -40.39068195549418, -40.391224188014434, -40.391766420534694,
                         -40.39230865305495,  -40.39285088557521,  -40.39339311809546,
                         -40.39393535061572,  -40.39447758313597,  -40.39501981565623,
                         -40.39556204817649],
                       dtype='float64', name='lat', length=2009))
        • lon
          PandasIndex
          PandasIndex(Float64Index([  -62.5120945327963, -62.511552300276044, -62.511010067755784,
                         -62.51046783523553,  -62.50992560271527,  -62.50938337019502,
                         -62.50884113767476, -62.508298905154504, -62.507756672634244,
                         -62.50721444011399,
                        ...
                         -61.24164373783563,  -61.24110150531537,  -61.24055927279512,
                         -61.24001704027486, -61.239474807754604, -61.238932575234344,
                         -61.23839034271409,  -61.23784811019384,  -61.23730587767358,
                        -61.236763645153324],
                       dtype='float64', name='lon', length=2353))
      • ncei_template_version :
        NCEI_NetCDF_Swath_Template_v2.0
        summary :
        The Earth Surface Mineral Dust Source Investigation (EMIT) is an Earth Ventures-Instrument (EVI-4) Mission that maps the surface mineralogy of arid dust source regions via imaging spectroscopy in the visible and short-wave infrared (VSWIR). Installed on the International Space Station (ISS), the EMIT instrument is a Dyson imaging spectrometer that uses contiguous spectroscopic measurements from 410 to 2450 nm to resolve absoprtion features of iron oxides, clays, sulfates, carbonates, and other dust-forming minerals. During its one-year mission, EMIT will observe the sunlit Earth's dust source regions that occur within +/-52° latitude and produce maps of the source regions that can be used to improve forecasts of the role of mineral dust in the radiative forcing (warming or cooling) of the atmosphere.\n\nThis file contains L2A estimated surface reflectances and geolocation data. Reflectance estimates are created using an Optimal Estimation technique - see ATBD for details. Reflectance values are reported as fractions (relative to 1).
        keywords :
        Imaging Spectroscopy, minerals, EMIT, dust, radiative forcing
        Conventions :
        CF-1.63
        sensor :
        EMIT (Earth Surface Mineral Dust Source Investigation)
        instrument :
        EMIT
        platform :
        ISS
        processing_version :
        V1.0
        institution :
        NASA Jet Propulsion Laboratory/California Institute of Technology
        license :
        https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/
        naming_authority :
        LPDAAC
        date_created :
        2022-11-14T09:50:54Z
        keywords_vocabulary :
        NASA Global Change Master Directory (GCMD) Science Keywords
        stdname_vocabulary :
        NetCDF Climate and Forecast (CF) Metadata Convention
        creator_name :
        Jet Propulsion Laboratory/California Institute of Technology
        creator_email :
        sarah.r.lundeen@jpl.nasa.gov
        creator_url :
        https://earth.jpl.nasa.gov/emit/
        project :
        Earth Surface Mineral Dust Source Investigation
        project_url :
        https://emit.jpl.nasa.gov/
        publisher_name :
        USGS LPDAAC
        publisher_url :
        https://lpdaac.usgs.gov
        publisher_email :
        lpdaac@usgs.gov
        identifier_product_doi_authority :
        https://doi.org
        flight_line :
        emit20220903t163129_o24611_s000
        time_coverage_start :
        2022-09-03T16:31:29+0000
        time_coverage_end :
        2022-09-03T16:31:41+0000
        software_build_version :
        010603
        product_version :
        01
        history :
        PGE Input files: radiance_file=/beegfs/store/emit/ops/data/acquisitions/20220903/emit20220903t163129/l1b/emit20220903t163129_o24611_s000_l1b_rdn_b0106_v01.img, pixel_locations_file=/beegfs/store/emit/ops/data/acquisitions/20220903/emit20220903t163129/l1b/emit20220903t163129_o24611_s000_l1b_loc_b0106_v01.img, observation_parameters_file=/beegfs/store/emit/ops/data/acquisitions/20220903/emit20220903t163129/l1b/emit20220903t163129_o24611_s000_l1b_obs_b0106_v01.img, surface_model_config=/beegfs/store/emit/ops/repos/emit-sds-l2a/surface/surface_20221020.json
        easternmost_longitude :
        -62.5120945327963
        northernmost_latitude :
        -39.3067591475017
        westernmost_longitude :
        -61.236221412633064
        southernmost_latitude :
        -40.39610428069674
        spatialResolution :
        0.000542232520256367
        spatial_ref :
        GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
        geotransform :
        [-6.25120945e+01 5.42232520e-04 -0.00000000e+00 -3.93067591e+01 -0.00000000e+00 -5.42232520e-04]
        day_night_flag :
        Day
        title :
        EMIT L2A Surface Reflectance 60 m V001

      Now we have an orthorectified dataset with coordinates that can be plotted. We can also mask the fill_values of -9999 to make clear figures.

      In [17]:
      out_xr = out_xr.where(out_xr['reflectance'] != fill_value)
      out_xr.isel(bands=b650).hvplot.image(cmap='viridis', aspect = 'equal', frame_width=500, rasterize=True)
      
      Out[17]:

      Contact Info:¶

      Email: LPDAAC@usgs.gov
      Voice: +1-866-573-3222
      Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹
      Website: https://lpdaac.usgs.gov/
      Date last modified: 01-04-2023

      ¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.